molmass_fructose = 180; % Molar mass of fructose
close all;

%% Do 25C fit
tempmat = monod_mat_fructose_25C;
smat = size(tempmat);

ivec = (1:4)';

conc1 = tempmat(ivec,1);
grmax1 = nanmean(tempmat(ivec,2:3),2);
grmax_err1 = nanstd(tempmat(ivec,2:3),0,2);

% Do fit on corrected concentration (in mM)
figure;
errorbar(conc1,grmax1,grmax_err1, 'bo', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'b', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);


%% Do 30C fit
tempmat = monod_mat_fructose_30C;
smat = size(tempmat);

ivec = (1:4)';

conc1 = tempmat(ivec,1);
grmax1 = nanmean(tempmat(ivec,2:3),2);
grmax_err1 = nanstd(tempmat(ivec,2:3),0,2);

% Do fit on corrected concentration (in mM)
errorbar(conc1,grmax1,grmax_err1, 'ko', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'k', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);

%% Do 37C
tempmat = monod_mat_fructose_37C;
smat = size(tempmat);

conc1 = tempmat(:,1);
grmax1 = nanmean(tempmat(:,2:3),2);
grmax_err1 = nanstd(tempmat(:,2:3),0,2);

% Do fit on corrected concentration (in mM)
errorbar(conc1,grmax1,grmax_err1, 'ro', 'linewidth', 1);
    
hold on;
smat = size(tempmat);

set(gcf, 'position', [0 0 400 300])
set(gca, 'fontsize', 20);
xlabel('Concentration (mM)')
ylabel('Growth rate (1/h)')

box off;

% Monod fit
conc_tot = [conc1];
grmax_tot = [grmax1];

monodfit = @(x,s) x(1).*s./(x(2) + s);
x0 = [1 0.01];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(monodfit, x0, conc_tot,grmax_tot);

hold on;
xfit = linspace(0,max(conc_tot), 1000);
plot(xfit, monodfit(x,xfit), 'r', 'linewidth', 1);

ci = nlparci(x,residual,'jacobian',jacobian);

err_x1 = abs(ci(1,2)-x(1));
err_x2 = abs(ci(2,2)-x(2));

disp('Fit and error of KM (mM)');
disp(x(2));
disp(err_x2);